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ABSTRACT 

We study, using both theory and shxiulations, a system of self-gravitating sheets. A 
new statistical mechanics theory — free of any adjustable parameters — is derived 
to quantitatively predict the final stationary state achieved by this system after the 
process of coUisionless relaxation is complete. The theory shows a very good agreement 
with the numerical simulations. The model sheds new light on the general mechanism 
of relaxation of self-gravitating systems and may help us to understand cold matter 
distribution in the universe. 
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1 INTRODUCTION 

One of the most fundamental questions of astrophysics con- 
cerns the mass distribution after a self-gravitating system 
has relaxed to equilibrium. This problem, posed more than 
6 decades ago, still remains unresolved, see the recen t papers 
l|Hiorth fc William j[2O10l : lBarnes fc Winiamj|201ll ) and ref- 
erences therein. In the astrophysical context, the answer to 
this question may help to shed light on a lot of intriguing the- 
oretical puzzles, such as the physical mechanism responsible 
for the regularities observed in the light profile of elliptical 
galaxies and the mass distribution in the dark matter halos. 

T he classica l simulations of iNavarro. Frenk fc Whit j 

l|l995l . ll996l . ll997l ) have produced the density profiles of dark 
matt er halos that no p resent theory is c a pable of explaining , 
see llCamml [l950,; Hohl fc FeiJ Il966l: iLvnden-Belll Il967j: 
Hohl fc Campbelll [ 19681: ICuperman. Goldstein fc LecaJ 
19691 : [Shul 1 19781 : IWhite fc NaravanI Il987l : ISaslawl I2OO0I ') 
and references therein. The scope of the problem extends 
from the foundations of statisti cal mechanics to the large 
scale evolution of the universe, (|Padmanabhanlll990l b The 
purpose of the present paper is to construct a statisti- 
cal theory that is able to successfully predict both the 
mass and the velocity distributions of a self-gravitating 
system after it has completed the process of coUisionless 
relaxation. In this respect the One Dimensional Self- 
Gravitating Model (ODSGM) of interacting mass sheets, 
has proven very useful for understanding both the stellar 
dynamics and t he cosmo l ogical mo dels , sec ( Wright, Miller 
fc Stein 1982; iMathuij Il99d : Ijovce fc WorrakitpoonponI 
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l2010al lbl: iMiller fc Roue^ I2OI0I : iJovce fc SicardI I2OIID 
and references therein. Although much simpler than 
the real three dimensional gravity, this model contains 
the essenc e of the gravitational proble m: long-range 
potential (|Campa. Dauxois fc Ruffol |2009| ) and collec- 
tive motion damped by particle-w ave interactions which 
leads to coUisionless relaxation. (iLevin. Pakter fc Teled 
2008al: iLevin. Pakter fc Rizzatd l2008bl : iTeles et al.l boTg 



Pakter fc LevinlbOllI ) and references therein. The ODSGM 



is particularly convenient because the sheet-sheet inter- 
action potential is n ot singular, so that the mode l can 
be easily simulated (|Noullez. Fanelli fc Aurelll 120031 '). At 
the same time it contains much of the same complicated 
statistical mechanics (Campa et al. 2009) of other systems 
with long-range forces: broken ergodicity and lack of mixing, 
which also plague gravity in three dimensions. 

Long before the empirical works of Navarro et al. 
(1995, 1996, 1997) Lynden-BeU (LB) proposed a statisti- 
cal theory of self-gravi tating systems based on the Vlasov 
equation jLvnden- Belli 11967). The LB theory has become 
known as the Theory of Violent Relaxation. One of the 
fundamental assumptions of LB was that violent relax- 
ation leads to an efficient phase-space mixing. The subse- 
quent simulations, however, have shown that this is not the 
case (Bindoni & Sccco 2007; Yamaguchi 20081: Levin et al. 
2008a.b: ITeles et all 120101 : INavarro et al.i boiOl '). Therefore, 
a new approach is needed if one wants to theoretically un- 
derstand the structure of self-gravitating systems in equilib- 
rium. 

We will work in thermodynamic limit — the number 
of sheets A'^ diverges, the mass of each sheet m goes to 
zero, while the total mass M of the system remains fixed, 
mN — M. The advantage of working with Id system is 
that the gravitational potential is unbounded from above 
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and bounded from below. As a consequence of this — af- 
ter a sufficient large time Tx , known as the Chandrasekhar 
time — a finite self-gravitating system must rela x to th e 
usual Maxwell-Boltzmann (MB) equilibrium (R vbickill97ll ). 
In the limit N ^ oo, the Chandrasekhar time diverges 
and the system becomes trapped in an out-of-equilibrium 



stationary state (IWright et alj 



Tsuchiva fc G ouda 



f 995', 'f 996' 



1982 



Reidl fc Miller 



2000: 



Yawn & Miller 



1995 



200a 



Teles et al.l |2010 ; Gu pta fc Mukamel 2010 ). Although 3d 



self-gravitating systems never relax to MB equilibrium they 
also become trapped in a quasi-Stationary State qSS (Levin 
et al. 2008a,b). In fact, for real 3d self-gravitating systems, 
such as elliptical galaxies, t he life time of the qSS is larger 
than their proper life span l|Binnev fc Tremainellioosl '). For 
ODSGM it has been observed that qSS state h as a very long 
life-sp an even for a small number of particles jWright et al.l 
[l98i. 

This letter is organized as follows. In section [5] we in- 
troduce the Id gravitational sheet model. In section |3] we 
comment briefly on its general dynamics and present the 
results of numerical simulation. In section 3] we derive the 
equation which is capable, over a short time span, to accu- 
rately describe the oscillations of the root-mean-square of 
particle positions. In section [S] we present the statistical- 
mechanics theory that is able to accurately predict both the 
velocity and the mass distribution in equilibrium, without 
any adjustable parameters. The more involved derivation of 
the envelope equation is left to the appendix of the paper, 

El 



2 THE MODEL'S DESCRIPTION 

The ODSGM consists of A'^ sheets of mass rrii in the y z 
plane, free to move along the a;-axis. To simplify the calcu- 
lations, we suppose that all the sheets have the same mass 
rrii = m = M/N, where AI is the total mass of the system. 
The sheets interact through the one dimension gravitational 
potential that satisfies the Poisson equation 



\7'^ip{x,t) = ATvGp{x,t) , 



(1) 



where G is the gravitational constant and p{x, t) is the mass 
density. 

We define dimensionless variables by scaling mass, 
lengths, velocities, the potential, the mass density and the 
energy with respect to A/, Lq (an arbitrary length scale), 
Vb = V27rGMLo, Vo = 2ttGMLo, po = M /2Lo and 
-Eo = MVq — 2ttGM^ Lo, respectively. All these scales cor- 
respond to setting G — M = 1 and to defining the dynamical 
time scale as, 

'1/2 



TD = (47rGpo) 



(2) 



It is then easy to show that a particle (sheet) at the ori- 
gin with a mass density p{x) — 5{x) creates a long-range 
gravitational potential. 



tl;{x) 



(3) 



Systems with long-range interaction are very peculiar be- 
cause in the thermodynamic limit the collision duration time 
diverges, and the dynamical evolution of the system is gov- 
erned exa£tl2_by_tlwcolIraon^ Boltzmann (Vlasov) equa- 
tion, (|Braun fc Hepdll977l ') 



Dt ot dx ov 



(4) 



where /(x,v,t) is the one particle distribution function, so 
that p(x, t) = J /(x, V, t)dv. Clearly real self-gravitating sys- 
tems do not rigorously obey the N ^ oo limit, nevertheless, 
in the astrophysical context the number of "particles" is 
usually very large, so that it constitutes a very good ap- 
proximation. 

The Vlasov equation has an infinite number of invari- 
ants called the Cassimirs (|Chavanij|200^ ). Two of these are 
mass and energy. 



dxdv /(x, V, t) = 1 
dxdv (^H-^ 



/(x, V,t) = So , 



(5) 



(6) 



respectively, where £o is the initial energy. Linear 
momentum is also a conserved quantity. However, by the 
symmetry of the problem and without loss of generality, 
it can be set to zero. In fact, any local functional of the 
distribution function is a Cassimir invariant of the Vlasov 
dynamics. 



3 VLASOV DYNAMICS 

Unlike coUisional systems which after a very short time re- 
lax to the Maxwell-Boltzmann distribution, the time evolu- 
tion of the Vlasov equation never ends. Instead it progresses 
to smaller and smaller length scales. In practice, however, 
both experiments and simulations have only a finite resolu- 
tion. It is in this, coarse-grained, sense that we say that the 
Vlasov dynamics evolves to the stationary state. Unlike the 
usual MB equilibrium, however, the stationary state reached 
through the process of coUisionless relaxation depends ex- 
plicitly on the initial particle distribution. 



3.1 Numerical Simulations 

Numerical solution of the Vlasov equation is very com- 
plicated. One alternative is to perform directly the 
N-body simulation in which all particles evolve over 
time according to the usual Newton equations of mo- 
tion (jNoullez. Fanelli fc Aiireiil 120031 ). We consider a 
ODSGM in which sheets are initially distributed in accor- 
dance with the water-bag distribution. 



fo{x,v) = l]Q{Xm - \x\)Q{Vm ~ \v\) 



(7) 



where O is the Heaviside step function and 77 = (4 XmVm)'^ 
is the normalization constant obtained from Eq. (O. The 
values Xm and Vm are the limits of the distribution function 
in the phase space, so that without loss of generality we 
set Xra = 1. Due to the symmetry of the distribution, the 
total linear momentum is null and the only quantity other 
than the total mass to be conserved is the total energy f = 
«™/6 + l/3. 

The halo formation is related to the initial condition 
through the only dimensionless parameter that remains after 
the rescaling of section [5] — the virial number, TL = 2T/U, 
where T is the total kinetic and U is the potential energy. 
When the virial condition 2T = [/ is not satisfied, TZ ^ 1, 
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Figure 1. The initial (a) and final (b) phase space of 10"' self- 
gravitating sheets initially distributed according to the water-bag 
distribution with the virial number TZq = 2.5. 



the imbalance between the kinetic and the potential energies 
causes the system to develop collective oscillations which 
are then damped by particle-wave interactions, leading to 
the core-halo phase separation, Fig.l. Our interest in study- 
ing the distributions with 7?. 7^ 1 arises from the fact that 
there is not even a qualitative theory that is able to ac- 
count for the mass and velocity distributio ns for such sys- 
tems lYamashiro. Gouda fc Sakagamilll992l '). 



4 THE TEST PARTICLE MODEL 

Systems with long-range forces display collective motion 
such as Langmuir waves in plasmas or Jeans waves in gravi- 
tational svstems. feinnev fc Tremaind (|2008l ). Since there are 
no collisions, interaction of individual particles with the den- 
sity waves is the mechanism that responsible for the transfer 
of energy. The particle-wave interaction dampens the den- 
sity waves while at the same time transfers large amounts of 
energy to the individual particles. The interaction is similar 
to a surfer "catching" a wave — some particles can gain a 
large amount of energy from the bulk oscillations to escape 
from the gravitational cluster and move to high energy re- 
gions of the phase space. These particles will then form a 
tenuous halo that will surround the dense low-energy core. 

To study the oscillation of the initial gravitational clus- 
ter, we define the envelope as Xe{t) = -^3 < x{ty >. Note 
that with this definition, at t = the envelope is exactly 
the same as the extent of the initial particle distribution 
Xe{0) = Xm- For sufhciently short times, most particles will 
still remain within the limits set by Xe{t). Differentiating 
Xe{t) twice with respect to time, we find an approximate 
dynamical equation satisfied by Xe{t), 



7^o 

Xe(t) 



(8) 



where TZo is the virial number of the initial distribution (the 
derivation is presente d in AppendixlX]). We shall call this the 
"enve lope equation" (|Wangler et al.lll99^ : lDavidson fc QuinI 
I2OOII ). Equation |5] enables us to study the dynamics of a 
test particle Xi, subject to the oscillating potential produced 
by a uniform density distribution delimited by —Xe{t) ^ 
X ^ Xs{t). A test particle moving in this potential will then 
evolve according to 



X^{t) = 



for \x,it)\ ^ x,{t) 
-sgn[xi{t) - Xe{t)] for \xi{t)\ ^ Xe{t) 



(9) 
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Figure 2. Test particle dynamics as compared to the A'-body 
dynamics simulation: (a) corresponds to the evolution of Xe{t) 
obtained from the A'^-body simulation (dots) while the solid line is 
calculated using the envelope equation. The time is in units of to- 

(b) Represents the phase space after the first few oscillations and 

(c) is the Poincare plot of the phase space produced by the test 
particles. Note the appearance of the resonance islands both in 
the N-body and test particle simulation. The test particle model 
allows us to accurately calculate the location of the resonant orbit 
with the maximum extent x^. This will determine the maximum 
energy attained by the halo particles. At t = 0, 15 test particles 
were uniformly distributed from the center to slightly beyond the 
maximum core radius Xm with initial velocity 0. 



where Xe{t) is given by Eq. ((8| and sgn is the sign function. 
In Fig. 2 we compare the complete A^-body simulation with 
the test particle trajectories obtained from the numerical 
solution of Eqs. ([8]) and The good agreement between 
the two reveals that for short times the envelope equation [8] 
is very accurate for describing the initial oscillations of the 
ODSGM. The formation of halo also occurs on a very short 
time scale — few oscillations of the envelope. Indeed, from 
Figs. 2 (b), (c), we see that the test particle model allows us 
to accurately locate the resonant orbits ()Gluckstern 1994 ). 
which also delimit the extent of the halo in the full A'^-body 
simulations. Thus, the test particle model allows us to calcu- 
late the maximum energy — corresponding to the resonant 
orbit — attained by any of the halo particles. This energy 
will be used in the next Section to construct the core-halo 
distribution function. 



5 THE CORE-HALO STATISTICAL THEORY 

The Jeans theorem states that "any steady-state solution of 
the coUisionless Boltzmann equation depends on the phase- 
space coordinates only through integrals of motion of the 
given potential, and any function of the integrals yields a 
steady -state solution of the coU isionless Boltzmann equa- 
tion" , (jBinnev fc Tremaindliool ') . Thus, any function of en- 
ergy is a solution to the Vlasov equation. In particular, the 
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Maxwell-Boltzmann distribution (MB) is also a stationary 
solution of the Vlasov equation. However, it is important to 
remember that MB is not a global attractor of the Vlasov 
dynamics — an arbitrary initial distribution of particles will 
not converge to the MB distribution, as happens for systems 
governed by short-range forces. In this paper we are inter- 
ested in calculating the distribution to which the system 
will evolve, starting from an arbitrary water-bag-like initial 
condition, after the process of coUisionless relaxation is com- 
pleted. 

We first observe that the Vlasov equation can be inter- 
preted as a convective derivative of density over the phase- 
space, so that the distribution function evolves as an incom- 
pressible fluid. In the case of initial water-bag distribution 
this means that the phase-space density cannot exceed that 
of /o, /(x, V, t) ^ 77. The mechanism of core-halo phase sepa- 
ration now becomes clear. Some particles enter in resonance 
with the rms bulk oscillations. They gain a lot of energy 
and escape from the main cluster, forming a tenuous halo. 
At the same time, evaporation of highly energetic particles 
dampens the core oscillations. This is similar to the pro- 
cess of evaporative cooling. The oscillations will stop when 
all of the free energy is exhausted. However, because of the 
constraint on the maximum density imposed by the Vlasov 
dynamics, the system cannot freeze — collapse to the mini- 
mum of the potential energy. Instead, the particle distribu- 
tion of the core approaches that of a fully degenerate Fermi 
gas. With this physical insight in mind, we now propose an 
ansatz solution to the Vlasov equation: 



fch{ 



= 779 [ep - e] + X© [eh - e] e [e - ep] 



(10) 



where ef is the Fermi energy, is the halo (resonance) en- 
ergy, and £{x,v) = + is the one particle energy. 
The halo energy is obtained simply by taking the most ex- 
tended position of test particle Xh when it crosses the v = 
axis, in others words en ~ \xh\- The Fermi energy, cf, delim- 
its the extent of the core region, and x measures the phase 
space density inside the halo. Both of these parameters are 
calculated self-consistently from the norm and energy con- 
servation equations, Eqs. ((5| and ([6]). 

Substituting Eq. [10] into the Poisson equation and inte- 
grating the core-halo distribution function over velocity, the 
dimensionless Poisson equation becomes 

{(V - x)V<^F - ip + xV^h - "0 for tp ^ep 
XV^^T^ for ^F^i^^th (11) 
for ?/) ^ eh . 

Solving this equation numerically and using the con- 
straints Eqs. ((U and (|6]) yield the complete solution for 
the one-particle distribution function. In Fig. 3 we compare 
the theoretically calculated density and velocity distribution 
functions to the full A''-body simulation. To obtain the dis- 
tributions in the simulations, the system is evolved for 1000 
dynamical times until the stationary state is achieved. To 
get good statistics the histograms are constructed by bin- 
ning the particle positions and velocities at four different 
times. An excellent agreement is found between the the- 
ory and the simulations, without any adjustable parameters. 
Note that the agreement is very good for both the virial 
numbers above and below TZo = 1. In the same figure we 
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Figure 3. Final position (a) and velocity (b) distributions for 
initial water-bag with TZo = 0.5 at t=1000 rp). (c) and (d) are the 
stationary distribution for the initial water-bag with TZq = 2.5. 
The dotted lines are the initial distributions. 



also plot, with the dashed lines, the initial density and the 
velocity distributions. It is clear that during the process of 
coUisionless relaxation the system moves far from the ini- 
tial condition. It is important to keep in mind that although 
the system has relaxed to the stationary state, this state 
does not correspond to the thermodynamic equilibrium. In 
particular, the velocity distribution does not have the equi- 
librium Maxwell-Boltzmann form. This also shows that one 
needs to take a special care when modeling coUisionless sys- 
tem using hydrodynamic-like equations, since there might 
not be thermodynamic equilibrium even locally. 



6 CONCLUSIONS 

We have presented a theory which allows us to a priori pre- 
dict the density and the velocity distribution functions for 
a system of gravitationally interacting sheets. The theory 
is quantitatively accurate without any fitting parameters. 
The theory and the simulations show that in the thermo- 
dynamic limit this system does not evolve to the usual MB 
equilibrium. Instead it becomes trapped in a non-ergodic, 
non-mixing state — once formed, the halo never again equi- 
librates with the core in the N ^ 00 limit. 

Although in this paper we have studied only water-bag- 
like initial distributions, the theory c an be easily gener alized 
to more complex functional forms (jTeles et al.l 12009 ). The 
next step in the development of the theory will be to include 
a specific cosmological model to account for the expansion 
of the universe. If this can be achieved, the theory might 
shed new light on the mass distribution in the dark-matter 
halos. 
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APPENDIX A: ENVELOPE EQUATION 



We define the "envelope" as = \/3 < >. Differentiat- 
ing twice with respect to time we find, 

3 < a::a3 > 3 < > 9 < a;i >^ 
Xe = h . (Al) 

AjQ J,e Jjg 

To simplify the first term, we suppose that the mass density 
oscillations are affine — preserve the uniform mass distribu- 
tion. In this case, 



dip 

< XX >— — < X—— > = 
dx 



Integrating, we find 

< XX 



1 n^'^ x\ 

-r— / dx. 

2Xeit) J_^^j^-i.-j Xe 



(A2) 



(A3) 



The second term of eq. (|A1I) involves only the mean 
values which, for short times, should remain close to those 
obtained using the initial particle distribution function /o. 



< XX >- 



< a; >= 



4 XmV„ 



a:: dx = 

2Vm I 3 

xdx / idx — 



(A4) 



Finally the envelope equation reduces to, 

Xe(t) = — — 1 



Xe{t) 



(A5) 



where TZo = -jj^ = Since Xe(0) = Xm = 1 if TZo = 1 then 
Xe{t) = 0, and the system does not develop macroscopic 
collective oscillations. 
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